{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def initialization(nP, dim, ub, lb):\n",
    "    Boundary_no = len(ub)\n",
    "\n",
    "    if Boundary_no == 1:\n",
    "        X = np.random.rand(nP, dim) * (ub - lb) + lb\n",
    "    else:\n",
    "        X = np.zeros((nP, dim))\n",
    "        for i in range(dim):\n",
    "            ub_i = ub[i]\n",
    "            lb_i = lb[i]\n",
    "            X[:, i] = np.random.rand(nP) * (ub_i - lb_i) + lb_i\n",
    "    return X\n",
    "\n",
    "def unifrnd(a, b, c, dim):\n",
    "    mu = (a + b) / 2\n",
    "    sig = (b - a) / 2\n",
    "    return mu + sig * (2 * np.random.rand(c, dim) - 1)\n",
    "\n",
    "def rndX(nP, i):\n",
    "    Qi = np.random.permutation(nP)\n",
    "    Qi = np.delete(Qi, np.where(Qi == i))\n",
    "    return Qi[:3]\n",
    "\n",
    "def runge_kutta(XB, XW, DelX):\n",
    "    dim = XB.shape[1]\n",
    "    C = np.random.randint(1, 3) * (1 - np.random.rand())\n",
    "    r1 = np.random.rand(dim)\n",
    "    r2 = np.random.rand(dim)\n",
    "\n",
    "    K1 = 0.5 * (np.random.rand() * XW - C * XB)\n",
    "    K2 = 0.5 * (np.random.rand() * (XW + r2 * K1 * DelX / 2) - (C * XB + r1 * K1 * DelX / 2))\n",
    "    K3 = 0.5 * (np.random.rand() * (XW + r2 * K2 * DelX / 2) - (C * XB + r1 * K2 * DelX / 2))\n",
    "    K4 = 0.5 * (np.random.rand() * (XW + r2 * K3 * DelX) - (C * XB + r1 * K3 * DelX))\n",
    "\n",
    "    XRK = (K1 + 2 * K2 + 2 * K3 + K4)\n",
    "    return 1/6 * XRK\n",
    "\n",
    "def run_optimizer(nP, MaxIt, lb, ub, dim, fobj):\n",
    "    Cost = np.zeros(nP)\n",
    "    X = initialization(nP, dim, ub, lb)\n",
    "    Convergence_curve = np.zeros(MaxIt)\n",
    "\n",
    "    for i in range(nP):\n",
    "        Cost[i] = fobj(X[i, :])\n",
    "\n",
    "    Best_Cost = np.min(Cost)\n",
    "    Best_X = X[np.argmin(Cost), :]\n",
    "\n",
    "    Convergence_curve[0] = Best_Cost\n",
    "\n",
    "    it = 1\n",
    "    while it < MaxIt:\n",
    "        it += 1\n",
    "        f = 20 * np.exp(-(12 * (it / MaxIt)))\n",
    "        Xavg = np.mean(X, axis=0)\n",
    "        SF = 2 * (0.5 - np.random.rand(nP)) * f\n",
    "\n",
    "        for i in range(nP):\n",
    "            _, ind_l = np.min(Cost), np.argmin(Cost)\n",
    "            lBest = X[ind_l, :]\n",
    "\n",
    "            A, B, C = rndX(nP, i)\n",
    "            _, ind1 = np.min(Cost[[A, B, C]]), np.argmin(Cost[[A, B, C]])\n",
    "\n",
    "            gamma = np.random.rand() * (X[i, :] - np.random.rand(dim) * (ub - lb)) * np.exp(-4 * it / MaxIt)\n",
    "            Stp = np.random.rand(dim) * ((Best_X - np.random.rand() * Xavg) + gamma)\n",
    "            DelX = 2 * np.random.rand(dim) * np.abs(Stp)\n",
    "\n",
    "            if Cost[i] < Cost[ind1]:\n",
    "                Xb = X[i, :]\n",
    "                Xw = X[ind1, :]\n",
    "            else:\n",
    "                Xb = X[ind1, :]\n",
    "                Xw = X[i, :]\n",
    "\n",
    "            SM = runge_kutta(Xb, Xw, DelX)\n",
    "\n",
    "            L = np.random.rand(dim) < 0.5\n",
    "            Xc = L * X[i, :] + (1 - L) * X[A, :]\n",
    "            Xm = L * Best_X + (1 - L) * lBest\n",
    "\n",
    "            vec = [1, -1]\n",
    "            flag = np.random.choice(vec, dim)\n",
    "            g = 2 * np.random.rand()\n",
    "            mu = 0.5 + 0.1 * np.random.randn(dim)\n",
    "\n",
    "            if np.random.rand() < 0.5:\n",
    "                Xnew = (Xc + flag * SF[i] * g * Xc) + SF[i] * SM + mu * (Xm - Xc)\n",
    "            else:\n",
    "                Xnew = (Xm + flag * SF[i] * g * Xm) + SF[i] * SM + mu * (X[A, :] - X[B, :])\n",
    "\n",
    "            Xnew = np.clip(Xnew, lb, ub)\n",
    "            CostNew = fobj(Xnew)\n",
    "\n",
    "            if CostNew < Cost[i]:\n",
    "                X[i, :] = Xnew\n",
    "                Cost[i] = CostNew\n",
    "\n",
    "            if np.random.rand() < 0.5:\n",
    "                EXP = np.exp(-5 * np.random.rand() * it / MaxIt)\n",
    "                r = int(np.random.uniform(-1, 2))\n",
    "\n",
    "                u = 2 * np.random.rand(dim)\n",
    "                w = unifrnd(0, 2, 1, dim) * EXP\n",
    "\n",
    "                A, B, C = rndX(nP, i)\n",
    "                Xavg = (X[A, :] + X[B, :] + X[C, :]) / 3\n",
    "\n",
    "                beta = np.random.rand(dim)\n",
    "                Xnew1 = beta * Best_X + (1 - beta) * Xavg\n",
    "\n",
    "                Xnew2 = np.zeros(dim)\n",
    "                for j in range(dim):\n",
    "                    if w[j] < 1:\n",
    "                        Xnew2[j] = Xnew1[j] + r * w[j] * np.abs((Xnew1[j] - Xavg[j]) + np.random.randn())\n",
    "                    else:\n",
    "                        Xnew2[j] = (Xnew1[j] - Xavg[j]) + r * w[j] * np.abs((u[j] * Xnew1[j] - Xavg[j]) + np.random.randn())\n",
    "\n",
    "                Xnew2 = np.clip(Xnew2, lb, ub)\n",
    "                CostNew = fobj(Xnew2)\n",
    "\n",
    "                if CostNew < Cost[i]:\n",
    "                    X[i, :] = Xnew2\n",
    "                    Cost[i] = CostNew\n",
    "                elif np.random.rand() < w[np.random.randint(dim)]:\n",
    "                    SM = runge_kutta(X[i, :], Xnew2, DelX)\n",
    "                    Xnew = (Xnew2 - np.random.rand(dim) * Xnew2) + SF[i] * (SM + (2 * np.random.rand(dim) * Best_X - Xnew2))\n",
    "\n",
    "                    Xnew = np.clip(Xnew, lb, ub)\n",
    "                    CostNew = fobj(Xnew)\n",
    "\n",
    "                    if CostNew < Cost[i]:\n",
    "                        X[i, :] = Xnew\n",
    "                        Cost[i] = CostNew\n",
    "\n",
    "            if Cost[i] < Best_Cost:\n",
    "                Best_X = X[i, :]\n",
    "                Best_Cost = Cost[i]\n",
    "\n",
    "        Convergence_curve[it-1] = Best_Cost\n",
    "        print(f\"Iteration {it}: Best Cost = {Convergence_curve[it-1]}\")\n",
    "\n",
    "    return Best_Cost, Best_X, Convergence_curve\n",
    "\n",
    "def sphere(x):\n",
    "    return sum(x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "# Example usage\n",
    "if __name__ == \"__main__\":\n",
    "    nP = 50\n",
    "    MaxIt = 500\n",
    "    lb = -100 * np.ones(30)\n",
    "    ub = 100 * np.ones(30)\n",
    "    dim = 30\n",
    "\n",
    "    Best_fitness, BestPositions, Convergence_curve = run_optimizer(nP, MaxIt, lb, ub, dim, sphere)\n",
    "\n",
    "    # Results\n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    plt.plot(Convergence_curve)\n",
    "    plt.title('Convergence curve')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Best fitness obtained so far')\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "csharp"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
